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A method for computing the thermopower in interacting systems is proposed. This approach, 
which relies on Monte Carlo simulations, is illustrated first for a diatomic chain of hard-point elas¬ 
tically colliding particles and then in the case of a one-dimensional gas with (screened) Coulomb 
interparticle interaction. Numerical simulations up to N > 10 4 particles confirm the general the¬ 
oretical arguments for momentum-conserving systems and show that the thermoelectric figure of 
merit increases linearly with the system size. 
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I. INTRODUCTION 

Thermoelectric materials are of interest due to their 
ability to convert waste heat into electricity by the See- 
beck effect or use electricity for cooling by the Peltier 
effect M. The efficiency of a thermoelectric material 
is a monotonously growing function of the dimensionless 
figure of merit, 


where T is the temperature, a is the electrical conduc¬ 
tivity, k is the thermal conductivity, and S is the ther¬ 
mopower (or Seebeck coefficient). Increasing ZT is a 
challenging task due to the interdependency of transport 
coefficients. In particular, the electrical conductivity and 
the electronic contribution to the thermal conductivity 
are related by the Wiedemann-Franz law Q, which, fol¬ 
lows at low enough temperatures from the single-particle 
Fermi liquid theory, states that the ratio aT /k is con¬ 
stant. Such limitations could be overcome, in principle, 
by the energy filtering mechanism , i.e., when trans¬ 
mission of electrons is possible only within a tiny energy 
window. 

On the other hand, very little is known about the ther¬ 
moelectric properties of interacting systems. In this case, 
analytical results are rare and numerical simulations face 
difficult problems. To numerically evaluate ZT , one may 
put the system into contact with two thermochemical 
baths (reservoirs) allowing for particle exchange with the 
system. The reservoirs are tuned at different temper¬ 
atures and electrochemical potentials in order to main¬ 
tain stationary particle and heat currents. By analyzing 
the response of these currents to the temperature and 
electrochemical potential difference, one can evaluate the 
transport coefficients and, in turn, ZT, with Eq. ©■ 
This method has been successfully applied to the one¬ 
dimensional (ID) dimerized gas of interacting hard-point 
particles M- In that model, due to the fact that particles 
only interact via instantaneous collisions, the particles in 


the reservoirs are in effect decoupled from those of the 
system. As such the reservoirs can be modeled as ideal 
gases so that the simulations can be facilitated greatly 
(see Ref. [Tlj for the detailed description of the algo¬ 
rithm). However, in more general systems with realistic 
interaction, the coupling between the reservoirs and the 
system is essential; the mere injection of particles from 
an ideal gas into the system may induce huge, unphysical 
interaction energy when an injected particle is too close 
to a system particle. 

Given this difficulty, which is unsolved yet to the best 
of our knowledge, we turn to the closed heat baths [3 El 
that only exchange heat with the system. We sandwich 
the system with two closed heat baths at different tem¬ 
peratures to establish a nonequilibrium setup. For inter¬ 
acting systems, the heat conductivity can be computed 
directly with this setup, but computing thermopower is 
challenging. In this paper, we solve this problem by the 
following steps: We first use the grand-canonical Monte 
Carlo method |l4| to compute the electrochemical poten¬ 
tial, p, as a function of the particle density, p , at a given 
temperature; then, with our nonequilibrium setup and 
by molecular dynamics, we compute the density differ¬ 
ence A p across the system, set in response to the tem¬ 
perature difference AT applied to the system. Finally, 
based on the established relation between p and p, we 
map the density difference into the thermoelectric volt¬ 
age difference, AV, so that the thermopower is computed 
as S = —AV/AT. As to the electrical conductivity a, 
due to the fact that closed heat baths do not support a 
charge current, it cannot be computed numerically with 
our nonequilibrium setup; However, it can be computed 
in equilibrium simulations by using the Green-Kubo for¬ 
mula fl5j . 

To test and illustrate the method for computing the 
thermopower in nonequlibrium simulations, we first con¬ 
sider the ID dimerized hard-point gas model [3]- Then, 
we numerically investigate the case of a ID gas of par¬ 
ticles with nearest-neighbor Coulomb interaction, mod¬ 
eling a screened Coulomb interaction between electrons. 
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Due to momentum conservation (more generally, due to 
the existence of a single relevant conserved quantity), we 
expect on general grounds 0 that the figure of merit 
ZT diverges in the thermodynamic limit, implying that 
the Carnot efficiency is reached in this limit. This result, 
so far illustrated by means of toy models 00 , is here 
confirmed in a more realistic model. Moreover, ZT ex¬ 
hibits a rapid, linear growth with the system size, which 
is a consequence of the recently reported Fourier-like be¬ 
havior of thermal conductivity [18l - |23| . 

II. NUMERICAL METHOD 

The equations connecting fluxes and thermodynamic 
forces within linear irreversible thermodynamics are [2~il . 
0 


jp 

ju 



L P u \ ( -V(^) \ 
L uu ) V V(3 )' 


( 2 ) 


where j p is the local particle current, j u is the local 
energy current, p is the electrochemical potential, and 
/3 = 1 /{ksT) is the inverse temperature (we set the 
Boltzmann constant k B = 1). The kinetic coefficients 
Lij (with i,j = p,u) are related to the familiar transport 
coefficients as 


_ „ r 1 detL c 

a= T Lpp ’ K= T^~ 0 7’ b = eT\L 

( 3 ) 

here e is the charge of each particle (set to be e = 1), and 
det L denotes the determinant of the (Onsager) matrix 
of kinetic coefficients. Thermodynamics imposes det L > 
0, L pp > 0, and L uu > 0, and the Onsager reciprocity 
relations ensure that L up = L pu . Following Eq. ©, the 
thermoelectric figure of merit thus reads 


J pu 




ZT = 


( L U p pLppY 

detL 


( 4 ) 


Hereafter we describe a method for the computation of 
the transport coefficients, and consequently ZT, in a 
generic interacting system. Note that, even though we 
consider the particle and energy flows along one direc¬ 
tion, the motion inside the system could be, in principle, 
two- or three-dimensional. 

Thermal conductivity.- We compute the thermal con¬ 
ductivity by nonequilibrium simulations. Two heat 
baths 0 0] at temperatures Tr = T + AT/2 and 
Tr = T - AT/2, respectively, are connected to the two 
ends of the system. Then the system is evolved and af¬ 
ter a long enough relaxation stage, when the stationary 
state has been established, the heat conductivity is evalu¬ 
ated as n = j u / (AT/ L), where the overbar denotes time 
averaging and L is the system size along the direction 
of the heat flow. The distributions of the temperature 
and the particle density are calculated at the stationary 
state as well. The temperature is numerically computed 


as T(x) = 2 ek(x)/p(x), where p(x) = 'j2id{x — Xi) is 

the particle density and ek(x) = JA S(x — x/) is the 
kinetic energy density &> . 

Thermopower.- We use the nonequilibrium setup with 
heat baths described above and prepare the system in 
the stationary state. The thermopower is defined as the 
magnitude of the induced thermoelectric voltage in re¬ 
sponse to the temperature difference across the system, 
i.e., S = —Ap/eAT, where Ap = p L — p R = eAV is the 
induced electrochemical potential difference. In our sim¬ 
ulations, we first compute the particle density pl and p R 
at the two ends of the system, and then map these val¬ 
ues into pl and p Rl respectively, by means of the grand- 
canonical Monte Carlo method (see the Appendix [A}. 

Electrical conductivity.- As we consider closed heat 
baths that do not exchange particles with the system, 
we cannot compute the electrical conductivity with our 
nonequilibrium setup. For this purpose, one can turn 
to equilibrium simulations by taking advantage of the 
Green-Kubo formula 0|. 


III. NUMERICAL TESTS WITH THE ID 
DIMERIZED GAS MODEL 

In order to test and elucidate the grand-canonical 
Monte Carlo method in computing the thermopower, 
we first consider a gas of N colliding hard-point par¬ 
ticles with alternate masses m and M, a paradigmatic 
model proposed in Re f. |_16 | and extensively investigated 
in the literatures 0, Il3l l23l. [H-0 for understanding 
low-dimensional transport problem. We start with the 
case of equal masses, m = M , for which the relation be¬ 
tween the electrochemical potential and the density is the 
same as for the one-dimensional ideal gas in the semiclas- 
sical limit [36j : 

p=k B Tln(pX). (5) 

Here, A = h/\/2nmk B T is the de Broglie thermal wave¬ 
length (h is the Planck’s constant). The excellent agree¬ 
ment between this analytical expression and the numer¬ 
ical Monte Carlo simulations is shown in Fig. [T] It 
is interesting to remark that from classical thermody¬ 
namics of a one-dimensional ideal gas we obtain p = 
k B T\w(Cp/\/T), where the constant C cannot be deter¬ 
mined by purely classical means. On the other hand, such 
ambiguity is not present in the grand-canonical Monte 
Carlo simulations. Indeed, this method is in some sense 
semi-classical, in that it uses as information the value of 
the de Broglie thermal wavelength and the grand canon¬ 
ical partition function where the particles are considered 
as indistinguishable (see the Appendix lAl the factor 1/A! 
in Eq. (IA1I) is of purely quantum origin). In our units, 
A = 1/VT, and therefore (7 = 1. 

As shown in Fig. [TJ with the grand-canonical Monte 
Carlo method we can numerically determine the depen¬ 
dence of the electrochemical potential on the particle den¬ 
sity at a given temperature T. Further, we can compute 
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FIG. 1: (Color online) Electrochemical potential versus par¬ 
ticle density for the equal-mass (m = M) hard-point gas. 
Numerical results of the grand-canonical Monte Carlo (MC) 
simulations (symbols) are compared with the analytical re¬ 
sults (lines) given by Eq. |5j. Here and in the following fig¬ 
ures we use units such that m = 1, the Boltzmann constant 
k _b = 1, and the de Broglie thermal wavelength A = 1/VT. 


the thermopower by nonequilibrium molecular simula¬ 
tions using two statistical thermal baths, with different 
temperatures Tl and T/j, coupled to the left and the right 
end of the system. When the first (last) particle collides 
with the left (right) side of the system, it is injected back 
with a new speed |u| determined by the distribution (37) 


Pl,r{v) 


IzbTl.r 


exp 


v 2 mi,N \ 

2fcsTi.fi/ ’ 


( 6 ) 


FIG. 2: (Color online) The temperature profile (top panel), 
the density profile (middle panel), and the thermopower (bot¬ 
tom panel) for the hard-point gas with unequal masses m = 1 
and M = (\/5-l-1)/2. The temperatures of the two heat baths 
are set to be Tl = 1.05 and Tr = 0.95. The dash-dotted line 
in the bottom panel shows the analytical result of S = 3/2. 


where mi and mjy are the masses of the first and the last 
particle. 

As an example, here we consider the gas of colliding 
hard-point particles with alternate masses m = 1 and 
M = (\/5 + l)/2 ss 1.618 [23j]. The mean distance be¬ 
tween two nearest-neighboring particles is set to be unity, 


Size 

TL 

pL 

Pl 

T'r 

PR 

Pr 

21 

1.033 

0.978 

-0.040 

0.979 

1.022 

0.032 

41 

1.037 

0.971 

-0.049 

0.967 

1.031 

0.046 

81 

1.041 

0.965 

-0.058 

0.960 

1.037 

0.054 

161 

1.045 

0.960 

-0.066 

0.956 

1.042 

0.061 

321 

1.047 

0.957 

-0.070 

0.953 

1.046 

0.066 

641 

1.049 

0.955 

-0.073 

0.952 

1.048 

0.068 

1281 

1.049 

0.954 

-0.074 

0.951 

1.050 

0.070 

2561 

1.050 

0.953 

-0.076 

0.951 

1.050 

0.070 

5121 

1.050 

0.953 

-0.076 

0.950 

1.051 

0.072 

10241 

1.050 

0.953 

-0.076 

0.950 

1.051 

0.072 


TABLE I: The numerically computed data for hard-point gas 
model. The temperatures of the two heat baths are set to be 
T l = 1.05 and Tr = 0.95. 


so that the system length (size) L equals the particle 
number N. Figure [2] shows the stationary temperature 
and density profiles (top and middle panels). The densi¬ 
ties pl and pr at the left and right ends of the chain can 
be computed directly. Then the relation between p and 
p provided by the grand-canonical Monte Carlo simula¬ 
tions allows us to obtain the corresponding values of pl 
and pr and to compute the thermopower (bottom panel 
of Fig. [2j as S = —Ap/eAT. It should be noted that, 
as shown in Table [0 for small system sizes, the internal 
temperatures T' L and T' R at the left and right ends of 
the system, at which the particle densities pl and pr are 
computed, are slightly different from the external tem¬ 
peratures Tl and Tr. This discontinuity is the result of 
a boundary resistance, generally denoted as Kapitza re¬ 
sistance; see for instance Ref. [l2| . This boundary effect 
vanishes when increasing the system size; see Table Q] 
However, it may be relevant for small systems sizes, as 
shown in the bottom panel of Fig. [2j where we compare 
the thermopower computed by using the internal tem¬ 
peratures as S = —(pl — pR)/e(T' L — T' R ) or the external 
temperatures as S = —(pl — Pr)/z{Tl — Tr). Since we 
are here interested in the intrinsic transport properties 
of the system rather than in the details of the coupling 
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to the heat baths, in what follows we will show results 
for the first method only. Moreover, we note that the 
first method shows a faster convergence to the asymp¬ 
totic analytical value S = 3/2 0 than the second one. 
Finally, the convergence of our numerical results for a 
nonintegrable model to the analytical value corroborates 
the validity of our computational scheme. 


IV. THERMOELECTRICITY OF THE 
COULOMB GAS 

With the help of the above described method, now we 
study a ID system of more general interaction. 


A. The model 



The model system we consider consists of N charged 
particles with nearest-neighbor Coulomb interaction, 
which is described by the Hamiltonian 


h = E 

i 




(7) 


where rrii, Xi, and pt = rniXi are the mass, the co¬ 
ordinate, and the momentum of the ith particle, re¬ 
spectively, and the nearest-neighbor interaction given by 
U(x) = a/x can be considered as a simplified effective 
model of screened interaction between particles (a is the 
controlling parameter of interaction strength and, due 
to Coulomb repulsion, particles cannot cross each other; 
i.e., Xi > Xi-\ and U(x) > 0). The overall momen¬ 
tum P = rriiXi is conserved. The total charge cur¬ 
rent is J e = eJ p , where J p = 'Yjii is the total particle 
current. The total energy current reads as follows fl2l ]: 
J u = Yhjii where the local energy current ji = ^(£i+i — 
Xi){x i+ 1 + Xi)F(x i+ 1 - Xi) + ±ihi with F(x ) = -U'(x) 
and hi = [rriiX~ + U{xi + \ — Xj) + U{xi — i)]/2. In 

the following we assume that all particles have the same, 
unitary mass; i.e., nrii = m = 1. Moreover, the mean dis¬ 
tance between two nearest neighboring particles is set to 
be unity, so that the system length L equals the particle 
number N. 


B. Simulation results 

We first compute by means of the grand-canonical 
Monte Carlo method the mapping between the density 
and the electrochemical potential for the Coulomb gas 
model. A broad value range of the interaction parameter 
a, ranging from a = 10 -4 to a = 1, has been investigated. 
Note that, as clearly shown in the top panel of Fig. [3j in 
the limit of a —> 0 the hard-point gas model is recovered. 
In addition, as expected and shown in the bottom panel 
of Fig. [3l the smaller a is, the closer the electrochemical 
potential (for a given particle density) to the analytical 


FIG. 3: Top: Coulomb potential U(x) = a/x for various 
values of the parameter a. Bottom: electrochemical potential 
versus particle density at temperature T = 1. In both panels, 
from top to bottom, a = 1,10 1 ,10 -2 , and ICC 4 , respectively. 
The straight dashed line shows the analytical relation [Eq. JoJ] 
between p and p for the hard-point gas. 


result given by Eq. © predicted for the hard-point gas 
model. 

For the computation of the thermopower and the ther¬ 
mal conductivity, we use Langevin heat baths [13] set at 
temperatures T B = T + AT/2 and T} j = T — AT/2. In 
our simulations, the value of AT, 10% of T, is set with 
the consideration that it is small enough to guarantee the 
system to be in the linear response regime but meanwhile 
not too small to facilitate the simulations. (Random tests 
with smaller AT, e.g., 4% of T have been done and the 
same results have been obtained.) The system is evolved 
with velocity-Verlet algorithm, but we have verified that 
all the results do not depend on the integration algo¬ 
rithm. The relaxation stage is longer than t = 10 7 for all 
the simulated cases. 

With regard to the electrical conductivity, in this case 
it is not necessary to use the Green-Kubo formula. In¬ 
deed, thanks to momentum conservation, a is ballistic 
and can be computed analytically. Using a stochas¬ 
tic model of thermochemical baths [H, [39j, we have 
j P = - 1r, where y a = p a y/k B T/\/2nm is the in¬ 

jection rate of particles from reservoir a into the sys¬ 
tem (a = L,R stands for the left and right reservoir, 
respectively), and p a is the density of particles in reser¬ 
voir a, modeled as an infinite one-dimensional ideal 
gas. Since the electrochemical potential p a for reser¬ 
voir a is given by p a = k B Thi(\p a ) [36|, we obtain 
7« = exp(/3/x a ) and, therefore, within linear re¬ 
sponse regime, j p = exp ^ M ° ; ) A p. The electrical con¬ 
ductivity is ballistic and given by er = GT, with the con- 
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x/L 



FIG. 4: (Color online) Temperature (top panel) and the den¬ 
sity (bottom panel) profiles for the Coulomb gas model. The 
temperatures of the two heat baths are set to be Tl = 1.05 
and T r = 0.95. 


ductance 


G= S = ^ eXp( ^ } ’ (8) 

where, always within linear response, /z ~ fiL ~ fJ-R- 
Figure |4] shows the stationary temperature and density 
profiles for the Coulomb gas model. Similarly to the di¬ 
atomic hard-point gas model, the system’s temperatures 
at the boundaries approach the temperatures of the ther¬ 
mal baths when the system size L — > oo; i.e., T' r — »• Tl 
and T' r —> Tr. In Fig. [5] we show the transport coef¬ 
ficients a, S, k, and the thermoelectric figure of merit 
ZT. For g we plot the analytically determined linear 
growth as a function of the system size, o = GL , with 
the conductance G given by Eq. © (we have checked 
that consistent results, i.e., g ~ L 1 can also be obtained 
by the Green-Kubo formula in equilibrium simulations, 
where the integration time is correctly truncated to take 
into account the ballistic transport [13 )■ The thermal 
conductivity increases for small system sizes and then 
saturates, as expected in a system that obeys the Fourier 
law. This behavior has been reported in recent investi¬ 
gations of several one-dimensional models of interacting 
particles ITsl l2dl . While the Fourier-like regime might be 
an intermediate (in the system size) regime, followed by 
an asymptotic regime of anomalous thermal conductiv¬ 
ity k ~ L 1 / 3 [12 . Il3j . its range may expand rapidly as 
an integrable limit (here, for a —> 0) is approached 12 .'ll . 
As a consequence, for practical purposes we can assume 
that the system obeys the Fourier law. Since also the 
thermopower saturates as predicted for ballistic trans¬ 
port [see Eq. (fl3l) in Sec. lIV C] , we can conclude that ZT 
grows linearly with the system size. In what follows, we 
show that the divergence of ZT with the system size can 




FIG. 5: (Color online) Dependence of the transport coeffi¬ 
cients a, S, k and ZT on the system size L for the Coulomb 
gas model. Dashed lines are drawn for reference. The tem¬ 
perature T — 1 and the interaction strength a = 1. 


be explained in terms of a general theoretical argument 
for momentum-conserving systems 0- 


C. Thermoelectricity in momentum-conserving 
systems 


At the thermodynamic limit, the presence of nonzero 
Drude weights T>ij is a signature of ballistic trans¬ 
port [40M43l ] : i.e., the kinetic coefficients Ljj scale lin¬ 
early with the system size L. As a consequence, the 
thermopower S is asymptotically size-independent. The 
finite-size Drude weights, for a system of size L , can be 
related to the existence of relevant conserved quantities 
of the system and computed by means of the Suzuki for¬ 
mula [44l l . Such formula states that 


Cij ( L ) = lim 


1 


t—> OO t 


dt' 


M 


= £ 


( JiQn)T(JjQn)T 

MYt 


( 9 ) 
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where (• • • )t denotes the thermal average at tempera¬ 
ture T, and {Q n , n = l,--- ,M} denote M orthogonal 
constants of motion, which are relevant; that is, non- 
orthogonal to the considered currents, in our case to the 
currents J p and J u : ( J P Q n )T 7 ^ 0 and ( J u Qu)t 7 ^ 0. The 
finite-size Drude weights are then defined as 

D ij {L) = -^jr Cij (L) . (10) 

If the thermodynamic limit L —> 00 commutes with the 
long-time limit t —> 00 , then the thermodynamic Drude 
weights T>ij can be obtained as 

T>ij = lim Dij(L). (11) 

L—too 

Moreover, if the limit does not vanish we can conclude 
that the presence of relevant conservation laws yields 
nonzero generalized Drude weights, which in turn implies 
ballistic transport. 

We can see from Suzuki’s formula that for systems with 
a single relevant constant of motion (M = 1), the bal¬ 
listic contribution to det L vanishes, since it is propor¬ 
tional to V pp V uu — V 2 , which is zero from Eqs. ©- 
m, and CD- Hence, detL grows slower than L 2 , and 
therefore the thermal conductivity n ~ detL/L pp grows 
sub-ballistically, n ~ L v , with v < 1. Furthermore, since 
a ~ Lgo ~ L is ballistic and S' ~ L°, we can conclude 
that jlq 

a Q 2 

ZT= - T oc L 1 ~ v . (12) 

K 

Hence, ZT diverges in the thermodynamic limit L —> 00 . 
This general theoretical argument applies, for instance, 
to systems where momentum is the only relevant con¬ 
served quantity. It has so far been illustrated in a toy 
model, i.e., a ID dimerized gas of interacting hard-point 
particles [lOj , and in a two-dimensional stochastic model 
of interacting particles E3- Here we consider the more 
realistic and complex model of the Coulomb gas. 

In order to check if our theory applies to the Coulomb 
gas model, in particular if the thermodynamic limit L —> 
00 commutes with the long-time limit t —> 00 , we need 
to compute the current correlation functions. For this 
aim we perform the equilibrium simulations with periodic 
boundary conditions. We prepare the equilibrium state 
of the system by using Andersen heat baths [45| at the 
same given temperature T and evolve the system for a 
sufficiently long time to make sure that it has been well 
thermalized. Then we remove the heat baths from the 
system. Starting from this moment the system is evolved 
isolatedly. After another long enough time of evolution, 
the correlation functions (<A(t) J, (0 ))t ( i,j = p,u) are 
computed as functions of time. 

The numerical results for the autocorrelation functions 
of J p and J U: and for the cross correlation function be¬ 
tween them are shown in Fig. [6] It can be seen that corre¬ 
lation function ( J u (t) J u {0))t approaches a finite nonzero 



t 

FIG. 6: (Color online) Current-current correlation functions 
for the Coulomb gas model at temperature T = 1 and inter¬ 
action strength a = 1. In all the panels, the black straight 
lines are for L = 64, the red dashed lines are for L = 128, 
and green dash-dotted lines are for L = 256. The results 
of the finite-size Drude weights Dij(L) by Suzuki’s formula 
(blue dotted horizontal line for L = 256) are also shown for 
comparison. 


value as the correlation time increases and that the char¬ 
acteristic time scale to approach such value is indepen¬ 
dent of the system size. Other current-current correla¬ 
tion functions are constant. We have, therefore, a strong 
numerical evidence that we can commute the long-time 
limit t —> 00 and the thermodynamic limit L —> 00 and 
we can compute the thermodynamic Drude weights 'D. lJ 
by means of Eq. CD- We also compute the finite-size 
Drude weights Dij(L) via the Suzuki formula Eqs. © 
and Eq. m, and indicate its value by a dotted hori¬ 
zontal line in the plots of Fig. [G] (Note that it is not a 
function of time t.) The obtained numerical results are in 
good agreement with the (asymptotical) values of the cor¬ 
relation functions. We note that ( J p {t) J p (0))r does not 
decay; this is because J p = P/m is a conserved quantity, 
so that the particle current is the same as for the ideal 
gas, and therefore by employing the Suzuki’s formula, 
we can analytically obtain that D pp = Tp/2. This result 
is in perfect agreement with the data shown in the last 
panel of Fig. El In fact, ( J p (t)J u (0))t = P{J u (0))t/ m 
is also trivially constant due to momentum conservation. 
(See the middle panel of Fig. [6]). Finally, as expected 
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from the theory 0, D pp (L)D uu (L) — D 1 2 pu (L) = 0; this 
is also verified for various system sizes. 

For momentum-conserving systems, we can also com¬ 
pute the asymptotic value of the thermopower (as L —> 
oo) based on theoretical prediction for ballistic transport: 


samples the grand-canonical probability distribution 
f ij.lt fe- N ] N ) oc L exp[-/?W(x JV )], (Al) 


S 


eT \ Dp P ) 


(13) 


Here, T>. t j can be obtained with the above equilibrium 
simulations via Eq. m, and c is a numerical con¬ 
stant that can be determined by comparison with the 
results obtained by the grand-canonical Monte Carlo 
method [46| . 


V. CONCLUSIONS 


where A is the de Broglie thermal wavelength and U is the 
potential energy for the iV-particle configuration x N = 
(xi, Xn). Our simulations are performed along the 
following steps (for a detailed description of the grand- 
canonical Monte Carlo method see Ref. 0): 

1. Start from an initial state with random positions of 
N particles; 

2. A random displacement is applied to a particle se¬ 
lected at random. This move is accepted with probability 


Starting from the definition of thermopower as a mea¬ 
sure of the magnitude of an induced thermoelectric volt¬ 
age in response to a temperature difference and taking 
advantage of the grand canonical Monte Carlo method 
to connect the particle density to the electrochemical po¬ 
tential, we are able to compute the thermoelectric coeffi¬ 
cients in systems with more general interaction than the 
instantaneous collisions. As a physically significant illus¬ 
tration of our approach, we have shown that for classi¬ 
cal one-dimensional, momentum-conserving systems with 
(screened) Coulomb interaction, the thermoelectric figure 
of merit increases, on a broad range of the system size, 
linearly. In principle, our strategy can be applied without 
restrictions on the type of interaction, even in the case 
of electron-lattice coupling, hence it could be useful for 
studying thermoelectricity in more complex and realistic 
systems. 


min{l, exp[-/I(Wnew - ^oid)]}, (A2) 

where W 0 id and U new denote, here and in the following, 
the potential energy before and after the move, respec¬ 
tively; 

3. The creation of a new particle at a random position 
is accepted with a probability 


^ + ^ exp[—/3(^ new — ^ 0 id)]|, (A3) 

where N a id denotes, here and in the following, the particle 
number before the move; 

4. The removal of a randomly selected particle is ac¬ 
cepted with a probability 
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Appendix A: Grand-canonical Monte Carlo 
simulations 

For a given electrochemical potential, system size, and 
temperature, the grand-canonical Monte Carlo method 


min < 1, 


5. Repeat steps 2 to 4, for a long enough time to reach 
the equilibrium state. 

6. Repeat steps 2 to 4, to have a sufficient number of 
microstates to compute the average number of particles 
(TV) and the density p = ( N)/L with good accuracy. 

Note that this algorithm obeys the detailed balance 
principle and therefore leads to a random sampling of 
the grand-canonical probability distribution 0. 


AAoid 


exp[ /3(W new ) 1 r > 


(A4) 
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